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1  Abstract 

A  realistic  simulator  of  complex  gas  phase  systems  needs  to  handle  accurately  processes 
spanning  different  scales  of  space  and  time.  The  paper  outlines  the  key  theoretical  and 
computational  features  of  modern  a  priori  treatments  of  the  dynamics  of  elementary 
processes  in  these  systems  with  particular  emphasis  on  reactive  processes. 


2  Introduction 

Realistic  a  priori  simulations  of  rarefied  gas  flows  are  conceptually  articulated  into  various 
blocks  of  operations  characterized  by  different  time  and  space  scales  (1).  These  blocks 
are  concerned  with  the  calculation  of  electronic  structures,  molecular  collisions,  collec¬ 
tive  (fluid-,  electro-,  etc.)  dynamics  and  the  averaging  over  randomly  sampled  variables 
to  work  out  a  priori  estimates  of  observed  and  measurable  properties  of  real  situations 
(like  the  formation  of  shock  waves  in  rarefied  gas  flows).  As  a  matter  of  fact,  these  are 
the  pillars  on  which  the  COST  Chemistry  Metalaboratory  SIMBEX  has  developed  the 
homonymous  simulator  of  molecular  beam  experiments  (1;  2)  and  the  COMPCHEM  (3) 
virtual  organization  (VO)  1  has  assembled  the  Grid  Enabled  Molecular  Simulator  (GEMS) 
proposed  for  European  funding  as  a  specific  targeted  research  project  within  the  Activity 
IST-2005-2.5.6  Research  Network  Testbed  of  the  6th  Framework  Program  (4). 

The  mentioned  blocks  of  computations  are  carried  out  using  separate  suites  of  codes 
thanks  to  the  fact  that  time  and  space  scales  associated  with  them  are  profoundly  dif¬ 
ferent.  This  means  that  the  electronic  energy  of  the  molecular  systems  and  the  related 
wavefunctions  can  be  calculated  to  a  high  level  of  accuracy  from  first  principles  using  ab 
initio  techniques  within  the  Born-Oppenhcimer  (BO)  (5)  approximation  (provided  that, 
when  necessary,  at  certain  times  of  the  process  or  at  certain  arrangements  of  the  nuclei  in 
the  corresponding  stationary  scheme,  the  coupling  between  nuclear  and  electronic  degrees 
of  freedom  is  regained).  In  the  BO  approximation  the  electronic  motion  is  assumed  to 
depend  only  parametrically  on  the  nuclear  coordinates  and  the  potential  energy  (elec¬ 
tronic  energy  plus  nuclear  repulsions)  of  the  nuclear  motion  is  calculated  (at  several  fixed 
geometries  of  the  nuclei)  using  well  consolidated  quantum  chemistry  suites  of  programs 
(see  for  example  refs.  (6;  7;  8;  9;  10))  which  will  not  be  discussed  here.  The  ab  initio 
calculation  of  the  electronic  energy  can  be  performed  either  time  by  time  at  the  actual 
geometry  of  the  moving  system  ( on-the-fly )  or,  once  for  ever,  at  the  first  step  of  the  com¬ 
putational  procedure.  In  the  latter  case,  ab  initio  calculations  are  performed  for  a  large 
matrix  of  nuclear  geometries  (including  those  of  the  initial  and  hnal  states  of  the  consid¬ 
ered  process)  and  calculated  values  are  then  best-fitted  by  optimizing  the  parameters  of 
a  suitable  potential  energy  functional. 

The  integration  of  the  motion  equations  of  the  nuclei  on  the  adopted  potential  energy 
surface  (PES)  allows  to  estimate  the  efficiency  parameters  of  the  elementary  chemical 
processes  being  considered.  These  detailed  efficiency  parameters  (with  their  temperature, 
energy,  spatial,  angular  momentum,  etc.  dependencies)  are  the  key  set  of  input  data 
needed  to  carry  out  realistic  simulations  of  the  kinetic  regime  of  real  gaseous  systems. 

In  this  paper  we  shall  discus  on: 

1COMPCHEM  is  registered  at  grid-it.cnaf.infn.it. 


RTO-EN-AVT-142 


10A-3 


Theory  and  Computing  of  Gas  Phase  Chemical  Reactions: 
From  Exact  Quantum  to  Approximate  Dynamical  Treatments 


ORGANIZATION 


1.  the  development  of  few  atom  rigorous  dynamical  methods  (section  2); 

2.  the  extension  of  dynamical  methods  to  larger  systems  (section  3); 

3.  the  exploitation  of  the  potentialities  of  advanced  computing  for  Molecular  Sciences 
calculations  (section  4). 


3  The  a  priori  dynamical  approach 


From  a  theoretical  point  of  view  elementary  chemical  processes  are  many  body  problems 
concerning  the  encounter  of  two  or  more  aggregates  of  electrons  and  nuclei.  This  allows  to 
treat  the  molecules  as  deformable  objects  which  can  collide  and  break  into  parts  and/or 
recombine.  Simultaneous  encounters  of  three  or  more  objects  are  quite  unlikely  to  oc¬ 
cur  (especially  in  the  low  pressure  gas  phase  processes  mentioned  before).  Accordingly, 
most  of  the  theoretical  and  computational  studies  have  focused  on  two  body  encounters 
(bi- molecular).  After  all,  the  treatment  of  multi- molecular  encounters  are  usually  ratio¬ 
nalized  by  chaining  relevant  bi-molecular  collisions  (though,  recently,  they  are  increasingly 
dealt  directly  using  classical  molecular  dynamics).  For  uni- molecular  processes  too  the 
theoretical  approach  can  be  easily  reconducted  to  that  of  half  a  bi-molecular  encounter. 

Typical  efficiency  parameters  of  elementary  bi-molecular  processes  are  either  (thermal, 
state-specific  and  state-to-state)  rate  coefficients  or  integral  (state-specific  and  state-to- 
state)  cross  sections.  The  thermal  rate  coefficient  k(T)  is  formulated  as  (11) 


*tn  =  EE 


Wj  exp  [- ej/kBT } 

Qint(T) 


h,f  (T) 


(1) 


where  i  and  /  are  respectively  the  initial  and  final  internal  states,  wy  is  the  multi¬ 
plicity  of  the  initial  internal  state  i  including  the  nuclear  spin  symmetry,  Qint(T)  = 
yh  Wi  exp  [— ei/ksT]  is  the  partition  function  associated  with  the  initial  internal  states, 
kg  is  the  Boltzmann  constant,  e*  is  the  energy  of  state  i,  kij(T)  is  the  state-to-state  rate 
coefficient  and  T  is  the  temperature  of  the  system.  The  evaluation  of  the  state-to-state 
rate  coefficient  can  be  reconducted  to  the  calculation  of  the  state-to-state  cross  section 
alj(Etr)  using  the  following  equation 


hf(T) 


EtrCJijiE^e-^/^dEtr 


(2) 


if  the  energy  distribution  is  of  the  Boltzmann  type  with  /j  being  the  reduced  mass  of  the 
system  in  its  reactant  arrangement  and  Etr  the  translational  energy.  The  state-to-state 
cross  section  ( ai}f(Etr ))  can  be  calculated  from  the  state-to-state  cumulative  reaction 
probability  (. Pij(Etr ))  as  follows 


d AEtr)  =  §P,:f(Etr)  =  §  E<2J  +  l)PiAE>r)  =  E<2' ’  +  !)  E  PU  (E‘r)  (3) 


J= 0 


J= 0 


A =— J 


where  ki  is  the  system  wavenumber  and  the  individual  terms  Pf^(Etr)  of  the  right  hand 
side  member  of  expression  3  can  be  derived  from  the  square  modulus  of  the  S  matrix 
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elements  calculated  at  a  given  pair  of  total  angular  momentum  quantum  number  J  and 
helicity  quantum  number  A. 


3.1  Electronically  adiabatic  potential  energy  surfaces 

As  already  mentioned,  a  priori  calculations  of  state-specific  and  state-to-state  probabil¬ 
ities  of  elementary  chemical  processes  start  from  the  consideration  that  the  fixed  nuclei 
electronic  wavefunctions  <f>e({w};  {W})  are  a  suitable  basis  set  for  the  expansion  of  the 
system  wavefunction  Z  ({w},  (W},f)  (in  our  notation  {w}  and  {W}  are  the  sets  of  elec¬ 
tronic  and  nuclear  position  vectors,  respectively,  and  t  is  time).  The  fixed  nuclei  electronic 
wavefunctions  <he({w};  {W})  are  the  eigensolutions  of  the  fixed  nuclei  electronic  problem 

He  ({w};  {W})  <f>e({w};  {W})  =  £({W})<T ({w};  {W})  (4) 

where  He  is  the  time  independent  electronic  Hamiltonian.  The  solution  of  eq.  4  provides 
us  with  the  desired  set  of  Jth  adiabatic  electronic  eigenvalues  £/({W})  for  the  molecular 
geometry  and  electronic  wavefunction  associated  with  it.  The  ensemble  of  the  £/({W}) 
values,  once  summed  to  the  corresponding  nuclear  repulsion,  represent  a  pointwise  de¬ 
scription  of  the  PES  (V/({W})  on  which  the  motion  of  the  Nnuci  nuclei  of  the  system 
takes  place. 

This  gives  us  a  means  for  solving  the  general  time  dependent  Schrodinger  equation  of 
the  system  (12;  13) 

({w},  {W},  t)  =  it {{ w},  {W })Z  ({w},  {W},  t)  (5) 

(in  eq.  5  H  =  He  +  Hn  is  the  total  many  body  Hamiltonian  of  the  system  whith  Hn  being 
its  nuclear  component)  by  expanding  Z,  as  already  mentioned,  in  terms  of  the  electronic 
eigenfunctions  which  parametrically  depend  on  the  position  vectors  of  the  nuclei  as  follows 


2  ({w},  {W},  t)  =  y  ¥,  ({W},  t)  $j  ({w};  {W}) .  (6) 

/ 

After  averaging  over  the  electronic  coordinates,  the  resulting  differential  equations 
for  the  T/({W},£)  coefficients  of  the  expansion  contain  some  terms  coupling  nuclear 
and  electronic  degrees  of  freedom.  As  already  mentioned,  the  BO  decoupling  scheme  is 
usually  applied  at  this  point  by  assuming  these  terms  to  be  negligible.  As  a  result,  for 
each  electronic  state  /  (hereafter,  the  index  /  will  be  dropped  from  the  formalism  because 
we  shall  confine  our  attention  to  the  single  electronically  adiabiatic  PES  BO  regime)  the 
calculations  reduce  to  the  following  electronically  adiabatic  time  dependent  Schrodinger 
equation  (12;  13;  14) 


ih^  ({W},  t)  =  hn({ W})T  ({W},  t)  =  [tw({w})  +  E({W})]  t  ({w},  t)  (7) 


where  T/v({W})  is  the  nuclear  kinetic  operator. 

In  regions  where  the  BO  approximation  breaks  down,  as  is  the  case  of  closely  spaced 
electronic  eigenvalues,  ad  hoc  treatments  can  be  adopted  which  deal  at  the  same  time 
with  different  terms  of  the  electronic  functions  manifolds  (15).  In  these  approaches  the 
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elements  coupling  electronic  and  nuclear  degrees  of  freedom  are  evaluated  consistently  us¬ 
ing  information  coming  from  electronic  structure  calculations.  Then  resulting  differential 
equations  coupling  different  electronic  states  and  wavefunctions  are  integrated  using  the 
same  standard  numerical  algorithms  available  to  integrate  equations  7. 

The  crucial  starting  phase  of  the  calculation  of  cross  sections  and  rate  coefficients  is 
therefore  the  assemblage  of  a  suitable  PES.  This  task  is  most  often  highly  demanding  in 
terms  of  computing  time  and  chemical  ingenuity.  For  this  reason  the  calculation  of  the 
potential  energy  values  is  increasingly  tackled  at  an  on-the-fly  ab  initio  level.  Yet,  up  to 
date,  the  most  popular  approaches  are  those  based  on  a  two  step  procedure.  The  first 
step  consists  of  the  collection  of  all  the  available  (both  ab  initio  and  empirical)  local  or 
global  data  on  the  system  interaction  while  the  second  step  consists  of  their  fitting  using 
an  appropriate  functional  form. 

For  small  systems  ( Nnuci  <  10),  as  is  the  case  of  the  majority  of  the  elementary  reac¬ 
tions  considered  up  to  date  for  theoretical  dynamical  studies  (in  some  cases  this  is  true 
also  for  larger  systems  if  the  complexity  of  the  computational  procedure  is  reduced  by 
imposing  suitable  dynamical  constraints),  the  above  mentioned  two  step  procedure  is, 
indeed,  the  preferred  one.  The  reason  for  this  preference  is  the  fact  that  the  quality  of 
present  electronic  structure  calculations  is  often  insufficient  to  guarantee  an  accurate  re¬ 
production  of  the  interaction  of  the  system  over  the  full  range  of  internuclear  distances 
unless  calculated  ab  initio  values  are  adjusted  using  empirical  considerations  before  un¬ 
dertaking  the  fitting.  The  functional  representation  to  be  used  for  fitting  the  PES  of 
reactive  systems  is  more  difficult  to  formulate  and  most  of  the  computations  have  been 
confined  to  systems  made  of  three  and  four  atoms  (16;  17;  18;  19;  20). 

The  most  popular  functional  forms  used  for  this  purpose  are  polynomials  either  in 
physical  coordinates  (17)  (like  the  internuclear  distances  defined  as  ry,-  =  |W,  —  Wj|)  or 
in  bond  order  variables  (21)  (defined  as  exponentials  of  the  displacement  from  equilibrium 
of  the  related  internuclear  distances).  When  using  physical  coordinates  the  polynomials 
need  to  be  damped  to  avoid  divergence  at  long  range.  Polynomial  functionals  are  usually 
adopted  within  a  Many  Body  Expansion  approach  to  formulate  the  individual  compo¬ 
nents  of  the  expansion.  Other  simple  functional  forms  are  either  derived  from  drastically 
simplified  formulations  of  ab  initio  methods  (22;  23)  or  from  intuitive  models  (such  as 
diatomic  rotating  potentials  (24;  25)).  Alternative  approaches  make  use  of  local  inter¬ 
polation  methods  in  which  for  each  interval  low  order  polynomials  are  employed  and  the 
value  of  related  parameters  are  determined  by  imposing  pointwise  or  switching  continuity 
conditions.  Similar  approaches  are  also  used  for  multi-surface  treatments  by  fitting  each 
surface  using  a  functional  form  (except  for  methods  directly  providing  multiple  solutions). 

When  moving  to  complex  systems  it  becomes  more  convenient  to  compose  the  PES 
by  summing  simple  few  body  functions  (stretches,  bends,  torsions,  van  der  Waals,  non- 
bonded  interactions,  etc)  containing  empirically  determined  parameters  (force  fields)  (26). 
Parameters  used  by  these  approaches  are  transferable  within  the  same  family  of  systems. 
Usually  these  surfaces  are  scarcely  suited  to  describe  the  making  or  breaking  of  bonds 
while  they  are  better  suited  for  conformational  studies. 
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3.2  Quantum  formalism  for  few  body  reactions 


Quantum  techniques  based  on  the  integration  of  eq.  (7)  are  widely  used  to  evaluate  the 
observable  properties  of  reactive  elementary  processes.  This  task  is,  nowadays,  a  largely 
routine  work  when  dealing  with  atom  diatom  systems.  From  the  integration  of  eq.  (7) 
one  can  evaluate  the  elements  of  the  scattering  matrix  S  whose  square  moduli  are  the 
elements  of  the  probability  matrix  P  used  in  eq.  3. 

For  a  few  body  isolated  system  (as  is  the  case  of  laser  free  crossed  molecular  beam 
experiments)  the  motion  of  the  center  of  mass  of  the  system  can  be  factored  out  to  reduce 
the  dimensionality  of  the  problem  (to  six  spatial  coordinates  in  the  atom  diatom  case 
that  is  the  simplest  prototype  of  elementary  reactions)  without  introducing  additional 
approximations. 

Accordingly,  for  a  three  atom  (A,  B,  C)  system,  the  time  dependent  Schrodinger 
equation  for  the  nuclei  reads  as 


d 

ih—  T(Rr,rr,f)  = 


,~(v'kr+vi)  +  v(Rr,rr,eT) 


T(Rr,r  T,t) 


(8) 


where  Rr  and  rT  are  the  mass  scaled  atom-diatom  Jacobi  vectors  (of  modulus  Rr  and 
ry),  ©T  is  the  angle  formed  by  Rr  and  rT  and  T(Rr,iy,f)  is  the  time  dependent  nuclear 
wavefunction.  As  usual,  Jacobi  coordinates  are  labeled  after  the  arrangement  r  (r  =  1,  2 
and  3  means  A  +  BC,  B  +  CA  and  C  +  BA  respectively)  to  which  they  refer. 

The  dimensionality  of  the  problem  is  further  reduced  if  the  Laplaeian  (the  kinetic 
component  of  the  Hamiltonian)  is  written  in  terms  of  angular  momentum  operators 


2/i 


(VL  +  vl) 


( J__^l ,  ,  A  ,  (J-Jt)2  _g_ 

2/i  \Rr  dR2  t  ry  dr2  T)  2 /iR2  2 fir2 


(9) 


with  J  being  the  total  angular  momentum  operator  given  by  the  sum  of  jT  and  lT  (the 
rotational  and  the  orbital  angular  momentum  operators  of  the  system,  respectively).  This 
makes  it  convenient  to  express  the  T(Rr,iy ,t)  wavefunction  in  terms  of  products  of  the 
xI>  JMp(R,T}  ry,  @T,  t)  partial  waves  (which  are  eigenfunctions  of  the  eigenvalue  J(  J+l)  of  the 
total  angular  momentum  operator  J2,  of  its  projection  M  on  a  space  fixed  (SF)  reference 
axis  and  of  the  parity  p)  in  the  internal  coordinates  RT1  ry  and  @r,  and  the  appropriate 
spherical  harmonics.  For  computational  convenience  to  integrate  the  scattering  equations 
from  the  reactant  to  the  product  arrangement  one  can  also  adopt  a  body  fixed  (BF) 
representation  (in  which  the  reference  frame  is  allowed  to  rotate  in  order  to  have  the  z 
axis  always  aligned  with  the  Rr  vector  and  the  xz  plane  having  a  fixed  orientation  with 
respect  to  the  molecular  plane)  where  A  is  the  projection  of  the  total  angular  momentum 
J  on  the  BF  z  axis. 


3.2.1  Time  independent  approaches 

To  further  reduce  the  dimensionality  of  the  problem  (this  was  the  approach  usually 
adopted  in  the  past)  the  time  variable  can  be  factored  out  from  the  system  wavefunc¬ 
tion  and  a  time  independent  formulation  of  the  Schrodinger  equation  can  be  obtained 
without  introducing  additional  approximations  (this  means  that  the  system  can  be  de¬ 
scribed  by  a  stationary  wave)  (27).  To  integrate  the  stationary  Schrodinger  equation  one 
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needs  to  define  a  particular  coordinate  (usually  called  reaction  coordinate)  by  properly 
combining  the  internuclear  distances  or  Jacobi  coordinates  to  ensure  a  smooth  switch 
from  reactant  oriented  to  product  oriented  arrangement  coordinates. 

The  most  popular  reaction  coordinate  adopted  in  recent  studies  has  been  the  hyper¬ 
radius  p  defined  as  p2  =  R2  +  r2  (this  relationship  holds  for  all  arrangements  because  p 
is  invariant  under  kinematic  rotations).  The  hyperradius  together  with  two  hyperangles 
constitutes  the  set  of  hyperspherical  coordinates.  These  coordinates  can  be  defined  in 
various  ways  depending  on  what  arrangement  is  to  be  preferred.  Here,  for  illustrative 
purposes,  we  make  use  of  their  APH  (28)  (democratic)  version  in  which  the  hyperangles 
are  9  and  x  (the  value  of  x  depends  on  the  choice  of  the  reference  geometry  though 
for  simplicity  the  related  label  is  dropped  here  from  the  notation)  and  the  partial  wave 
equations  take  the  form 


Tp  +  Th  +  Tr  +  Tc  +  V  (p,  9,  x) 


<tJM’’(p,e,x)  =  E<aJM’'(p,e,x ) 


(10) 


where  subscripts  “h”,  “r”  and  “c”  stand  for  “hypersphere”,  “rotational,”  and  “Coriolis”, 
respectively  and  the  operators  Tp,  Th ,  Tr  and  Tc  are  formulated  as: 


T  = 

1p 


h2  d  5d_ 
2/ip5  dpP  dp' 


h2  (  4  d  .  d  1  d2  \ 
“  2pp^  \skT26d6Sm  20 d9  +  sin2  9  dX2 J 

Tr  =  A(p,  9)J2  +  B(p,  9)J2  +  C(p,  9)Jl 


and 


-  ih  cos  9  d 
pp2sm29  v  dx 

with  A{p19)1  B(p,9)  and  C(p,9 )  being  defined  as  A^l{p,9)  =  /ip2(  1  +  sin  0),  B^1{p19)  = 
2pp2  sin2  9 ,  C-1(p,  9)  =  pp2(  1  —  sin  9). 

Eq.  (10)  is  integrated  by  segmenting  the  hyperradius  in  several  sectors  and  expanding 
T  (within  each  sector  i)  in  terms  of  the  surface  functions  which  are  eigensolutions  of 
the  following  equation 


15  h2 
8 


+  {C-D )  k2 A2  +  Dh2J{J  +  1)  +  V(Pi,  9,  X) 


(11) 

with  D  =  (A  +  B)/2  (though  the  $  functions  could  be  chosen  also  to  be  independent  of 
J  by  setting  J  =  0  in  eq.  11  and  regaining  the  J  dependence  during  the  integration  over 
p).  Once  the  expansion  is  performed  one  gets  the  following  set  of  equations  to  integrate 
over  the  hyperradius  p 


~  d2 
dp2 


2p 

¥ 


t'  A' 


where  the  internal  Hamiltonian  Hi  is  defined  as 


(12) 


Hi  —  B  +  /,•  +  Tc  + 


15  h2 


V(p,9,X). 


(13) 
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Accordingly,  the  computational  procedure  can  be  partitioned  into  three  parts.  Part 
one  is  devoted  to  the  calculation  of  the  $  functions  over  suitable  fixed  p  grid  points  of  9  and 
X  values.  This  part  is  dominated  by  the  evaluation  of  (two  dimensional  for  three  atom,  five 
dimensional  for  four  atom  systems)  integral  quadratures  and  eigenvalues  finding  of  large 
real  dense  square  matrices.  The  second  part  of  the  calculation  consists  in  propagating 
the  ^ln(Pi)  solutions.  This  part  of  the  calculation  is  dominated  by  the  inversion  of  large 
matrices  having  the  same  dimension  as  the  number  of  channels.  The  third  part  is  less 
computationally  demanding  and  is  devoted  to  the  mapping  of  the  asymptotic  solution 
into  the  proper  arrangement  space  and  to  the  imposition  of  boundary  conditions  to  work 
out  the  S  matrix. 

Moving  from  three  to  four  atom  systems  (say  diatom-diatom)  the  number  of  internal 
coordinates  doubles.  A  possible  choice  is  the  set  of  coordinates  A,  7q,  r2,  @1,  ©2  and  0 
with  R  being  the  distance  between  the  centers  of  mass  of  the  two  diatoms,  rq  and  r2  being 
the  two  diatomic  internuclear  distances,  ©i  and  ©2  being  the  two  planar  angles  formed 
by  rq  and  r2  with  R  and  0  being  the  dihedral  angle  formed  by  the  ( R ,  rq)  and  (R,  r2) 
planes.  Using  these  coordinates  the  Hamiltonian  of  the  AB  +  CD  system  takes  the  form 
(29) 


Hn  = 


h2  d2 
2/i4  dR2 


+  hi(rq)  +  h2(r2)  + 


+ 


jI 


+ 


j'i 


2 /i4i?2  '  2pi r\  '  2 /r2r, 


+  AD 


(14) 


where  /i4  is  the  reduced  mass  of  the  AB  and  CD  reduced  masses  (pi  and  /i2,  respectively), 
J  is,  as  already  mentioned,  the  total  angular  momentum  operator,  j42  is  the  sum  of  j\ 
and  j-2  which  are  the  rotational  angular  momentum  operators  of  AB  and  CD,  respectively, 
h\{ri)  and  h2(r2)  are  the  vibrational  Hamiltonians  of  AB  and  CD,  respectively,  while  A 
V  is  the  difference  between  the  total  interaction  potential  V(R,  rq,  r2,  ©i,  ©2,  0)  and  the 
V(rq)  and  U(r2)  vibrational  potentials  of  /ii(rq)  and  /r2(r2). 

The  increased  complexity  of  the  four  atom  Hamiltonian  makes  the  definition  of  the 
reaction  coordinate,  the  calculation  of  the  sector  basis  functins  and  the  switch  from  one 
arrangement  to  another  (and  therefore  the  solution  of  the  time  independent  Schrodinger 
equation)  very  difficult.  As  a  matter  of  fact,  only  recently  significant  advances  have  been 
made  in  describing  the  reaction  coordinate  of  four  atom  systems  by  using  row-orthonormal 
hyperspherical  democratic  coordinates  made  of  a  hyperradius  and  five  hyperangles  (30). 


3.2.2  Time  dependent  approaches 

The  difficulty  of  handling  the  problem  associated  with  the  definition  of  a  suitable  smoothly 
evolving  spatial  continuity  variable  brings  the  discussion  back  to  the  decision  of  reduc¬ 
ing  the  dimensionality  of  the  calculation  by  factoring  out  the  time  dependence  of  the 
wavefunction.  In  this  respect  the  simplicity  obtained  in  formulating  the  computational 
machinery  when  using  the  Jacobi  coordinate  time  dependent  formalism  is  increasingly 
considered  a  suitable  reward  for  keeping  the  extra  variable  of  time  in  the  formalism.  As  a 
matter  of  fact,  in  time  dependent  approaches  one  has  the  advantage  of  straightforwardly 
shaping  the  initial  wavepacket.  This  is,  in  fact,  chosen  to  correspond  to  a  given  reac¬ 
tant  state  or  to  a  mixture  of  them,  its  component  along  the  atom  diatom  coordinate  R 
is  formulated  as  a  product  of  a  gaussian  wave  (exp  [—a(R  —  -R0)2])  times  a  phase  shift 
factor  (exp  [— ik0(R  —  i?D)])  associated  with  its  inward  traveling  nature  having  an  average 
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momentum  kQh.  Then  the  wavepacket  can  be  mapped  into  any  other  set  of  arrangement 
coordinates  of  interest  and  let  propagate  in  time  by  repeatedly  applying  the  evolution 


operator  exp 


— iHr/h 


Eventually,  after  the  wavepacket  has  spread  all  over  the  whole 


accessible  configuration  space  one  can  carry  out  its  analysis  at  the  product  asymptotic 
line  by  expanding  the  cut  of  the  wavefunction  into  the  related  diatomic  wavefunctions 
(for  atom-diatom  systems).  The  time  dependent  coefficients  of  the  expansion  Cv'f\>(t ) 
read  (31) 


Cv'j'A'{t)  —  j  J  sin  @PJA(0)0„/J/(r)'kJA(i?  =  r,  0,  t)d0dr  (15) 

when  using  the  reactant  Jacobi  coordinates  R ,  r  and  0  (primed  variables  are  for  products, 
unprimed  for  reactants).  In  eq.  (15)  PA(0)  is  the  normalized  associated  Legendre  function 
of  the  angular  part  of  the  wavefunction.  By  performing  a  half  Fourier  transform  of  the 
coefficients  one  gets  the  time  independent  (energy  dependent)  A  matrix  whose 
elements  read 

AvtjiA'{E)  =  —  /  exp  (iEt/h)  ■  CvijiA'(t)dt.  (16) 

2vr  Jt=0 

From  them  one  can  easily  determine  the  S  matrix  elements  whose  square  modulus  (the 
probabilities)  allow  to  calculate  the  atom-diatom  cross  section  (31) 


(Tvj.v'j'  (Ejr  ) 


k2  (2 7  +  11  ^  \^vjA,v'j'A' 

Kvj  +  )  j  p=±l  A, A' 


(17) 


3.2.3  Direct  calculations  of  Reactions  rate  coefficients 

The  simplicity  of  the  time  dependent  method  has  also  facilitated  its  extension  to  larger 
systems  (29;  32).  However,  when  one  is  interested  in  evaluating  the  (less  detailed)  thermal 
rate  coefficient  there  is  no  need  to  carry  out  in  the  computation  all  the  details  embodied 
in  the  S-matrix.  After  all,  the  value  of  the  rate  coefficient  is  rather  insensitive  to  the 
detailed  structure  of  the  whole  PES  while  it  is  strongly  dependent  on  the  shape  of  the 
reactant  side  of  the  saddle  and  in  particular  on  the  height  and  the  width  of  the  reaction 
barrier  (33). 

This  has  motivated  the  formulation  of  the  rate  coefficient  directly  in  terms  of  the 
cumulative  reaction  probability  N(E)  defined  as 

=  (18) 

i  f 

and  that  does  not  refer  to  any  asymptotic  state  and  depends  only  on  the  dynamics  of  the 
system  in  the  vicinity  of  the  reaction  barrier.  In  terms  of  N(E)  the  rate  coefficient  is  then 
formulated  as 

1  r°° 

k{T)  =  l  N^-EltBTAE  <19> 

where  Qtr(T )  is  equal  to  (27 T/iksT)3^  /h3.  This  means  that  the  efficiency  of  the  reactive 
process  is  expressed  in  terms  of  the  fraction  of  the  system  wavefunction  left  over  after 


10A- 10 


RTO-EN-AVT-142 


Theory  and  Computing  of  Gas  Phase  Chemical  Reactions: 
From  Exact  Quantum  to  Approximate  Dynamical  Treatments 


projecting  out  all  its  components  which  do  not  have  outgoing  character  in  the  product 
asymptotic  region.  The  projection  operator  can  be  given  the  form 

Pprod  =  lim  Fm  6{R  -  R0)  .  (20) 

t—>  OO 


Accordingly,  the  wavefunction  can  be  propagated  forward  (out  of  an  arbitrarily  located 
dividing  surface)  infinitely  in  time.  Then  components  not  located  on  the  product  side  of 
the  dividing  surface  are  projected  out  and  the  remaining  wavefunction  is  propagated 
backward  in  time.  In  this  approach,  the  rate  constant  k(T )  can  be  formulated  as: 


k(T) 


1 

Qtr(T)Qint{T) 


lim  tr(F  e-A/kBT  elftt  h  e~iAt) 

t— >oo 


(21) 


when  such  a  limit  exists.  In  equation  21,  h  can  be  chosen  to  be  any  operator  discriminating 
between  reactants  and  products  and  F  any  operator  measuring  the  flux  from  reactants 
to  products  (in  Eq.  (21)  the  following  correlation  function 

Cfp(t )  =  tr(e~^r  F  e~^r  elflt  h,  e~lAt)  (22) 


of  the  flux-position  type  is  used).  As  an  alternative  use  can  be  made  of  the  following 
correlation  function 

Cff(t)  =  tr(e~ air  fe'®'  elflt  F  e~lAt)  (23) 

of  the  flux-flux  type  with  Cff(t )  being  the  time  derivative  of  Cfp(t )  and  F  the  time 
derivative  of  h  in  the  Heisenberg  picture.  Accordingly,  the  thermal  rate  coefficient  takes 
the  form 


k(T)  = 


lim  Cfp(t)  = 


Qtr(T)Qint(T )  ^  '  Qtr(T)Qint(T)  J 

o 

and  the  cumulative  reaction  probability  becomes  (34;  35): 

N(E)  =  2t r2  tr  (f  5{H  -  E)  F  5{H  -  E) 

1 


cffm 


(24) 


2C 


2  E/kBT 


St  S't 


AEt 


<  fr 


- iHt 


fT  >  dt 


(25) 


where  the  evaluation  of  the  trace  (whose  detailed  calculation  would  imply  the  propagation 
of  the  whole  set  of  basis  functions)  has  been  obtained  in  terms  of  the  eigenstates  fr  of 
the  thermal  flux  operator  Ft  defined  as 

Ft  =  e~A/2kBT  F  e~A/2kBT  (26) 

which  need  only  a  small  number  of  applications  of  Ft  on  a  sample  wavefunction. 

An  interesting  analogy  between  the  classical  and  the  quantum  description  of  the  rate 
coefficient  can  be  obtained  if  one  factorizes  the  correlation  factor  into  a  static  and  a 
dynamic  component.  The  dynamic  component  provides,  in  fact,  in  the  t  — >  oo  limit, 
a  description  of  the  amount  of  Pt  eigenstates  ending  up  on  the  product  side  of  the 
dividing  surface.  Different  ways  of  exploiting  these  ideas  to  the  end  of  carrying  out  the 
actual  calculations  of  the  value  of  the  thermal  rate  coefficients  and  of  cumulative  reaction 
probabilities,  are  given  in  refs.  (36;  37;  38;  39). 
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4  Approximate  Methods 

The  practical  impossibility  of  carrying  out  exact  calculations  for  more  complex  systems 
has  prompted  the  development  of  computational  procedures  based  on  various  decoupling 
schemes  effective  in  reducing  the  dimensionality  of  the  problem. 

4.1  Reduced  dimensionality  quantum  methods 

The  goal  of  working  out  approximate  formulations  of  the  cross  section  and  of  the  rate 
coefficient  was  first  accomplished  for  atom  diatom  systems  by  introducing  decoupling 
schemes  of  the  kind  energy  sudden  (40),  centrifugal  sudden  (41)  and  infinite  order  sudden 
(42).  In  these  decoupling  schemes  either  diatomic  rotations  or  atom  diatom  orbiting  (or 
both)  are  treated  in  an  approximate  way.  For  systems  made  of  four  or  more  atoms  the 
coupling  of  the  various  degrees  of  freedom,  while  significantly  increasing  the  complexity 
of  the  dynamical  treatment,  most  often  plays  a  negligible  role  in  determining  the  reaction 
outcome. 

This  has  allowed  a  split  of  the  dynamical  treatment  of  strongly  coupled  degrees  of 
freedom  from  weakly  coupled  ones.  Strongly  interacting  degrees  of  freedom  are  treated 
rigorously  while  weakly  interacting  ones  are  treated  approximately.  The  most  popular 
approximations  are  based  either  on  adiabatic  assumptions  or  on  the  parametrization  of 
some  variables.  As  an  example,  in  the  rotating  bond  approximation  (RBA)  of  diatom- 
diatom  non  linear  collisions,  two  radial  coordinates  and  one  bending  angle  are  explicitly 
treated  while  the  other  three  degrees  of  freedom  are  kept  frozen  (43).  The  overall  result 
is  then  worked  out  by  averaging  quantities  calculated  at  the  different  values  of  the  frozen 
variables.  In  the  adiabatic  bend  approximation  (ABA)  (44;  45),  instead,  the  three  radial 
coordinates  are  treated  explicitly  while  the  three  bending  angles  are  treated  adiabatically. 
This  means  that  the  overall  wavefunction  is  factored  out  and  eigenvalues  associated  with 
the  effect  of  the  terms  of  the  Hamiltonian  in  the  related  coordinates  are  used  to  obtain 
the  effective  Hamiltonian  for  the  coordinates  to  be  treated  exactly. 

As  already  mentioned  the  calculation  of  dynamical  properties  can  be  performed  by 
determining  the  cumulative  probability  which  can  be  estimated  in  an  approximate  way 
using  a  transition  state  (TS)  schematization  of  the  reactive  process.  In  this  view  the  role 
played  by  the  overall  rotation  of  the  system  is  that  of  shifting  in  energy  (of  a  quantity 
E shift  related  to  the  rotational  energy  of  the  system  at  the  TS  geometry)  the  reaction 
probability  P-Jj(Etr).  This  approximation,  usually  called  J-shifting  (45),  links  the  value  of 
the  state-to-state  probability  calculated  at  a  given  value  of  the  total  angular  momentum 
quantum  number  J  to  that  calculated  at  a  reference  value  ( Jref  )•  The  relationship  used 
is 

Pijj,  (E,r)  =  F$.,  (Etr  -  E 3hi ft)  (27) 

where  Eshift  is  the  difference  in  energy  between  the  overall  rotational  eigenstates  J  and 
Jref  at  the  transition  state  and  Jref  is  usually  taken  equal  to  zero.  A  more  elaborated 
method  linearly  interpolates  P/j(Etr )  between  the  probabilities  calculated  at  the  two  Jref 
values  most  closely  sandwiching  the  actual  J  value  of  the  calculation. 
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4.2  Mixed  quantum  and  classical  mechanics  approaches 

Approadies  alternative  to  quantum  calculations  are  those  based  on  classical  mechanics. 
Several  of  these  approaches  use  classical  mechanics  formulations  to  describe  slower  mo¬ 
tions  and  quantum  mechanics  formulations  to  describe  the  remaining  (faster)  degrees 
of  freedom.  Related  schemes  have  been  applied  to  divide  electronic  from  nuclear  mo¬ 
tion  and  to  formulate  widely  general  scaling  procedures  (47;  48).  Some  of  these  mixed 
quantum-classical  approximations  have  been  derived  in  a  rigorous  way  in  refs.  (49;  50)  by 
introducing  a  specific  basis-set  of  the  orthogonal  polynomial  type  centered  at  a  ”  classical” 
trajectory  and  expanded  around  it  as  the  dynamics  evolves  in  time.  In  this  case  the  forces 
appearing  in  the  classical  equations  of  motion  are  not  the  usual  Newtonian  ones  but  more 
general  forces  usually  called  ”  quantum  forces” . 

The  chosen  basis  can  be  used  in  combination  with  the  usual  time  independent  basis 
functions  or  collocated  representations.  This  allows  to  work  out  mixed  quantum  classical 
approaches  and  to  monitor  the  quantum  classical  correlation  and  measure  the  accuracy 
of  classical  path  treatments. 

When  systems  become  very  large,  the  computational  procedures  need  to  be  simplified 
further  also  at  the  level  of  the  calculation  of  the  electronic  energies  (in  this  case,  in  fact,  the 
accuracy  of  the  calculated  electronic  wave  function  becomes  intrinsically  poor  (51)).  As  a 
result,  the  choice  of  election  is  to  combine  the  use  of  classical  mechanics  for  treating  the 
dynamics  of  the  nuclei  with  the  adoption  of  approximate  schemes  to  calculate  electronic 
energy  (like  in  density  functional  approaches  (51))  and  greatly  simplify  the  on-the-fly  (52) 
computational  machinery. 

These  Ab-initio  Molecular  Dynamics  simulations  are  presently  (53;  54)  applied  to 
the  study  of  physico-chemical  properties,  such  as  reactivity  and  dynamic  relaxation,  of 
several  systems  (55).  Hybrid  procedures  linking  Quantum- Mechanical  parameterized  de¬ 
scriptions  of  the  ‘active  sites’  of  the  molecule  with  a  Molecular  Mechanics  description  of  its 
inactive  framework  are  highly  popular.  These  Quantum  Mechanics  Molecular  Mechanics 
(QMMM)  treatments  (56)  show  advantages  and  limitations.  They  still  lead,  in  fact,  to 
dynamical  results  which  cannot  be  inferred  by  a  mere  static  analysis  of  the  features  of 
the  potential  energy  surface  though  they  do  not  treat  most  of  them  accurately. 

At  various  stages  of  a  trajectory  calculation  one  has  the  possibility  of  treating  semiclas- 
sically  those  degrees  of  freedom  for  which  a  classical  approach  is  inadequate  by  associating 
to  them  a  semiclassical  wave  depending  on  the  classical  action  accumulated  along  the  clas¬ 
sical  path.  This  allows  to  regain  concepts  like  flux,  interference,  resonance  and  tunneling 
within  a  trajectory  framework  and  reproduce  quite  closely  some  quantum  features  of  the 
results. 


4.3  The  pure  classical  mechanics  computational  machinery 

Purely  classical  mechanics  approaches  found  the  dynamical  treatment  on  the  integration 
of  classical  trajectories  (CT).  For  this  reason  they  are  called  quasiclassical  or  QCT  when 
initial  and  final  states  are  in  some  way  discretized.  CT  methods  assume  that  the  nuclei 
involved  in  a  chemical  reaction  obey  classical  mechanics  and  roll  as  point  mass  particles  on 
the  potential  energy  surface  of  the  system.  Accordingly,  H^t(TI),  the  classical  analogue 
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of  the  quantum  nuclear  Hamiltonian  Hn( R)  of  eq.  (7),  reads 

Nnucl  rrp  —I—  Tp  —I—  Tp 

R})  =  £  ^  Pk^+V({R}) 


(28) 


and  the  evolution  of  the  system  is  followed  by  integrating  the  equations  of  classical  me¬ 
chanics  starting  from  different  sets  of  initial  conditions  (for  atom  diatom  systems  initial 
conditions  are  given  by  the  vibrational  number  is,  the  rotational  number  j,  the  transla¬ 
tional  energy  Etr,  the  impact  parameter  b,  the  elongation  of  the  diatom,  the  phase  of 
the  rotating  diatom,  the  angles  formed  by  the  rotational  angular  momentum  and  the 
velocity  with  the  molecular  plane).  Various  formulations  of  the  classical  particles  motion 
equations  can  be  given.  In  the  widely  used  Hamilton’s  version  they  read  as 

ds^  _  d Hn 
d t  dpSRk 

&PsRk  =  dHN 
d  t  dsRk 

for  each  cartesian  projection  sRk  (of  the  position  vector  R*.)  and  pSR  (of  the  momentum 
vector  pk)  of  each  atom  k  of  mass  my.  of  the  molecular  system.  The  equations  are  inte¬ 
grated  starting  from  one  set  of  allowed  initial  conditions  of  the  reactants  in  state  i  and 
are  stopped  either  when  the  maximum  number  of  interactions  steps  has  been  reached  or 
when  an  asymptotic  geometry  of  the  products  has  been  reached.  As  already  mentioned, 
discrete  features  of  quantum  results  are  then  enforced  in  QCT  ones  by  arbitrarily  dis¬ 
cretizing  the  energy  of  classical  bound  motions.  Usually  for  atom  diatom  systems  this 
means  that  the  counter  NVjUij>  (associated  with  trajectories  starting  from  isj  and  ending 
with  a  classical  vibrotational  energy  closer  to  that  of  the  ls'j'  state  than  to  any  other  one) 
is  incremented  by  a  unit.  When  all  planned  Ntj  trajectories  are  integrated  PUjyj'  is  set 
equal  to  Nvjyy  /Ntj  and  < 7Vjyy  is  set  equal  to  nbmaxPUjyj'  (with  bmax  being  the  maximum 
value  used  for  the  impact  parameter).  CT  and  QCT  approaches  can  often  provide  esti¬ 
mates  of  rate  coefficients,  cross  sections,  angular  distributions  and  reaction  probabilities 
of  reasonable  accuracy.  Moreover,  they  allow  a  pictorial  view  of  the  mechanisms  governing 
chemical  reactions.  Obviously,  the  CT  method  is  an  approximation  to  the  nuclear  motion 
and  it  becomes  more  accurate  when  quantum  effects  are  negligible  (as  is  the  case  of  heavy 
nuclei,  large  collision  energies  and  highly  averaged  reactive  properties). 


5  Advanced  Computing 

As  it  has  been  already  mentioned,  the  demand  for  computer  resources  prompted  by  chem¬ 
ical  reactivity  calculations  and  related  realistic  a  priori  simulations  is  as  high  as  that  of 
other  grand  challenges  of  modern  computational  science.  To  guarantee  the  computational 
feasibility  of  these  applications  one  has  to  resort  to  the  exploitation  of  the  innovative 
features  of  parallel  and  distributed  computers  by  decomposing  the  problem  in  simpler 
decoupled  subproblems  and  distributing  the  resulting  independent  blocks  of  the  codes 
for  execution  on  a  large  quantity  of  processing  elements.  The  application  of  some  de¬ 
compositions  based  on  physical  considerations  (like  separating  the  electronic  structure 
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calculation,  factorizing  the  time  variable,  disentangling  the  center  of  mass  motion  and 
performing  a  partial  wave  expansion)  have  been  already  discussed  before.  However,  when 
dealing  with  large  systems  and  complex  applications  one  has  to  effectively  exploit  the  in¬ 
novative  features  of  parallel  and  distributed  computing  by  carrying  out  a  decomposition 
of  the  problem  at  an  algorithmic  level. 

5.1  Parallel  computing 

To  effectively  tackle  the  problem  of  parallel  restructuring  a  computational  application  the 
following  aspects  need  to  be  specifically  considered  (57;  58;  59) 

1.  the  key  algorithms  need  to  be  optimized  for  parallel  execution, 

2.  existing  software  modules  need  to  be  integrated  efficiently, 

3.  different  programming  models  and  languages  need  to  be  used  at  various  levels, 

4.  retrieving  and  storing  of  data  structures  need  to  be  reorganized  and  streamed, 

5.  performances  of  the  adopted  articulation  need  to  be  measured  under  different  con¬ 
ditions  for  its  improvement. 

This  particularly  difficult  and  time  consuming  job  is  usually  carried  out  by  reorganizing 
the  relevant  suites  of  codes  and  by  inserting  the  appropriate  directives  and  commands 
chosen  among  those  of  the  most  popular  parallelization  libraries  (60).  However,  the  need 
for  ensuring  both  reliability  and  standardization  on  one  side  and  the  difficulty  of  keeping 
the  pace  of  the  continuous  evolution  of  architectures  and  simulation  techniques  on  the 
other  side  have  made  it  necessary  to  produce  tools  guaranteeing  the  automatic  or  semi¬ 
automatic  portability  of  applications  (also  in  the  sense  of  performance  portability)  across 
computing  platforms. 

Significant  progress  along  this  direction  has  been  made  using  structured  environments. 
A  typical  structured  environment  useful  for  a  semi-automatic  parallelization  of  the  appli¬ 
cations  is  SKIE  (61).  SKIE  is  an  integrated  environment  providing  a  new  application  ori¬ 
ented  set  of  instruments  allowing  the  rapid  development  and  prototyping  of  applications. 
Such  an  environment  is  based  on  some  optimized  and  ready-to-use  parallel  structures, 
called  skeletons.  The  skeletons  can  embody  sections  of  codes  allowing  so  far  an  extended 
reuse  of  the  existing  sequential  (written  in  the  most  popular  high  level  sequential  lan¬ 
guages)  or  parallel  programs  by  encapsulating  them  in  modules.  Examples  of  skeletons 
are  processor  farms  (a  pool  of  worker  processes  computes  a  pool  of  independent  tasks); 
pipelines ,  (different  processes  carry  out  in  a  sequence  the  various  phases  of  the  computa¬ 
tion);  map  data-parallel  computations  (all  the  elements  of  a  data  structure  are  updated 
or  computed  at  the  same  time).  These  structures  are  handled  using  a  coordination  lan¬ 
guage  (CL)  called  SKIE-CL  and  can  be  utilized  to  coordinate  and  connect  any  sequential 
or  parallel  module  encapsulated  using  a  SKIE-CL  wrapper.  The  wrapper  ensures  that 
parameter  passing  and  data  representation  are  consistent  among  the  modules  composing 
a  parallel  application.  SKIE-CL  makes  use  of  instruments  like  control ,  stream  parallel 
and  data  parallel  though  it  accepts  also  in  the  input  and  output  parameter  lists  all  the 
usual  basic  types  of  variables  (integers,  real,  etc),  records,  and  multi-dimensional  arrays. 
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In  addition,  SKIE-CL  provides  a  stream  data  type  which  allows  programmers  to  generate 
or  consume  sequences  of  items  of  indefinite  length.  The  peculiar  feature  of  these  patterns 
is  the  fact  that  they  can  be  freely  composed  ”a  la  carte”  to  build  complex  structures. 

SKIE  automatically  generates  also  an  optimized  implementation  of  the  skeleton  com¬ 
position.  This  means  that  when  using  a  SKIE  skeleton  the  support  not  only  generates 
automatically  the  code  needed  for  parallel  interaction,  but  it  also  optimizes  the  resources 
allocated  to  each  skeleton,  decides  the  best  granularity  of  computation  and  locates  ineffi¬ 
ciencies  in  the  global  structure. 

Going  to  a  more  abstract  level  skeletons  can  be  generalized  practically  to  any  form 
(or  combination  of  forms)  and  size.  This  is  indeed  the  key  feature  of  the  coordination 
language  ASSIST  (62;  63)  that  is  an  evolution  of  SKIE.  ASSIST  is  made  of  a  graph 
whose  nodes  are  parallel  ( parrnod )  and  sequential  ( seqmod, )  modules.  The  arches  of  the 
graph  represent  streams  of  data.  Interaction  among  parmods  can  occur  also  via  shared 
objects  (data,  memory  regions,  functions,  etc.)  for  which  ASSIST  provides  a  primitive 
access  mechanism  by  implementing  or  emulating  a  shared  memory  access.  Non  primitive 
accesses  are  instead  provided  for  external  objects  (like  DSM  libraries,  CORBA  servers, 
etc.)  for  which  the  user  has  to  take  care  of  the  access,  synchronyzation  and  consistency 
via  directives  to  be  inserted  in  the  code.  Activities  in  different  parmods  can  be  parallel 
or  concurrent  (like,  for  example,  in  a  pipeline).  Parallel  or  concurrent  activities  can  take 
place  also  within  the  same  parrnod.  They  can  be  either  farm-like  or  data-parallcl-likc  (or 
mixed  in  a  non- deterministic  way  depending  on  the  structure  of  the  data  and  the  status  of 
the  computation).  This  means  that  ASSIST  allows  the  following  two  hierarchical  levels: 
among  various  parmods  and  within  the  same  parrnod.  The  first  level  describes  a  graph 
of  the  data  flow  type  while  the  second  can  describe  computations  both  of  a  data  parallel 
domain  decomposition  type  and  of  a  functional  replication  farm  type.  ASSIST  is  also 
scarcely  invasive  since  the  computation  is  described  using  a  set  of  procedures  wrapping 
the  sequential  user  code  and  organized  in  the  parallel  fashion  specified  by  the  used  parrnod 
type. 

ASSIST  has  been  already  used  to  parallelize  time  dependent  reactive  scattering  ap¬ 
plications.  Applications  implemented  using  ASSIST  and  its  libraries  (see  ref.  (64))  have 
shown  to  clearly  outperform  their  versions  implemented  using  MPI.  Similar  studies  have 
been  carried  out  for  time  independent  applications  separately  for  the  surface  functions 
calculations  and  the  propagation  along  the  reaction  coordinates.  Tests  performed  using 
MPI  on  an  eigenvalue  finding  routine  are  illustrated  in  ref.  (65) 

5.2  Virtual  organizations  and  Grid  enabled  applications 

Moving  to  realistic  simulations  of  complex  chemical  systems,  however,  parallelism  on  a 
single  (no  matter  how  powerful)  machine  is  still  insufficient.  As  a  matter  of  fact,  the 
emerging  computing  paradigm  is  the  computing  Grid  (66;  67).  The  Grid  offers  the  pos¬ 
sibility  of  solving  complex  problems  using  a  (preferably  very  large)  set  of  distributed 
computers  as  a  single  unified  computing  resource.  To  this  end  the  Grids  enable  the  shar¬ 
ing,  selection,  and  aggregation  of  a  wide  variety  of  geographically  dispersed  resources 
ranging  from  PCs  to  supercomputers,  storage  systems,  data  sources  and  specialized  de¬ 
vices.  These  may  be  owned  by  different  organizations  and  work  for  completely  different 
purposes  within  the  virtual  organization.  Accordingly  a  Grid  can  be  viewed  as  a  seamless 
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integrated  computational  platform. 

Moreover,  the  Grid  is  by  definition  at  the  same  time  both  a  collaborative  environment 
and  a  market  in  which  the  users  interact  to  find  a  solution  for  their  problems  and/or  to  offer 
their  solutions  to  other  users.  Just  to  mention  a  few  services  which  can  be  ’’traded”  on  the 
Grid  we  can  list:  Computational  services,  Data  services,  Application  services,  Information 
services  and  Knowledge  services.  This  makes  Grid  applications  intrinsically  large  scale 
and  multidisciplinary.  Obviously,  in  order  to  make  the  Grid  computing  effective,  a  large 
number  of  tools  concerning  security,  resource  allocation,  costs  management,  information 
flux,  software  development,  process  execution,  resource  aggregation  and  scheduling  need 
to  be  implemented.  As  a  matter  of  fact  the  Grid  is  able  to  use  in  a  synergistic  way 
resources  of  all  kinds  including  those  which  cannot  be  physically  replicated  in  a  single 
site.  This  means  that  one  can  easily  scale  up  computing  cycles  as  well  as  competences 
regardless  of  their  location  to  work  in  a  coordinated  fashion.  In  summary,  the  Grid  is  a 
virtual  place  for  composition  of  multiple  administrative  domains  and  autonomies  to  work 
concertedly  on  a  large  variety  of  heteregenous  machines  and  problems  which  naturally 
scales  up  from  a  few  integrated  resources  to  very  many  of  them  a  smooth,  dynamic, 
adaptable  and  interoperable  way. 

The  European  Grid  of  election  is  the  already  mentioned  EGEE  (68).  Within  EGEE 
a  memorandum  of  understanding  has  been  signed  last  March  with  the  COMPCHEM  vir¬ 
tual  organization  to  foster  the  development  of  molecular  science  complex  simulations.  The 
mapping  of  a  complex  application  on  the  computing  grid  is,  indeed,  a  difficult  task.  Such 
a  process  is  not  simple  because  there  is  not  a  unique  correspondence  between  the  variation 
of  physical  and  mathematical  parameters  of  a  complex  computational  application  and  the 
distribution  on  the  grid  of  the  computing  blocks.  This  is  not  only  due  to  the  fact  that 
a  variation  of  the  computational  parameters  alters  the  relative  importance  of  the  various 
computational  blocks  and  of  the  related  demand  of  computational  resources  but  also  to 
the  fact  that  the  support  of  the  grid  infrastructure  at  run  time  is  not  deterministic  (69). 
For  this  reason  it  is  vital  to  work  out  data  graphs  and  build  workflow  managers  allowing 
a  proper  independent  handling  of  the  computing  blocks  at  various  levels  of  distribution. 
The  Erst  step  of  this  process  is,  therefore,  the  breaking  of  the  computational  procedure 
into  independent  or  loosely  coupled  computing  blocks.  In  fact,  the  singling  out  of  inde¬ 
pendent  (or  almost  independent)  computational  tasks  is  propedeutic  to  the  design  of  any 
concurrent  organization  of  the  relevant  computer  programs. 

5.3  The  Grid  implementation  of  a  molecular  science  simulator 

The  particular  Grid  enabled  application  considered  as  a  study  case  by  COMPCHEM  is, 
as  it  has  been  already  mentioned,  GEMS  the  Grid  Enabled  Molecular  Simulator.  A  demo 
version  of  GEMS  (GEMS.O)  has  been  already  implemented  on  the  production  Grid  of 
EGEE  and  presented  at  the  first  EGEE  review  workshop  (70).  GEMS.O  is  derived  from 
the  SIMBEX  simulator  (see  ref.  (1))  developed  by  the  homonimous  working  group  of  the 
COST  Chemistry  Action  D23  (Metalaboratories  for  complex  computational  applications 
in  Chemistry)  (71).  It  takes  care  of  evaluating  the  cross  sections  and  the  product  distri¬ 
butions  (plus  some  non  observable  quantities)  of  a  crossed  molecular  beam  atom-diatom 
experiment  using  a  quasiclassical  approach.  In  a  quasiclassical  approach,  the  observable 
properties  of  a  scattering  experiment  are  determined  by  performing  a  multidimensional 
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integral  over  the  unselected  variables  of  the  experiment.  The  kernel  of  the  integral  is 
given  by  a  boolean  function  whose  value  depends  on  the  result  of  the  integration  of  the 
Hamilton  equations  of  the  molecular  system. 

This  means  that  one  has  to  integrate  large  batches  (say  N tj)  of  classical  trajectories 
for  which  some  of  the  initial  conditions  are  randomly  selected.  The  calculation  of  each 
individual  trajectory  implies  the  integration  of  the  equations  of  motion  of  the  Nnuci  atoms 
composing  the  system.  This  means  that  the  simulation  has  to  generate,  possibly  in  a 
deterministic  way,  a  large  number  of  subsets  of  pseudorandom  numbers  (with  each  subset 
determining  the  initial  conditions  of  a  given  trajectory).  The  overall  workflow  of  GEMS.O 
is,  therefore,  articulated  into: 

1.  a  first  part  defining  the  parameters  of  the  calculation,  computing  quantities  of  gen¬ 
eral  use  and  assembling  the  information  needed  for  the  calculation  of  the  potential, 

2.  a  final  part  performing  the  averaging  of  the  calculated  quantities  and  carrying  out 
the  graphical  elaboration  of  the  properties  to  be  rendered  possibly  in  real  time  on 
the  virtual  monitors, 

3.  a  central  (the  key)  part  iterating  over  the  distribution  to  the  workers  of  the  trajec¬ 
tories  to  be  integrated  and  to  the  recollection  of  the  related  results  to  update  the 
reactive  probability. 

The  central  part  of  GEMS.O  is  the  kernel  of  the  calculation  that  can  be  efficiently  dis¬ 
tributed  on  the  Grid  using  a  task-farm  scheme.  The  distribution  can  take  place  by  as¬ 
signing  the  integration  of  individual  (or  blocks  of)  trajectories  to  the  worker  nodes  after 
generating  initialseed  (the  first  seed  of  the  subset  needed  by  each  trajectory  or  subset 
of  trajectories)  at  Master  process  level.  In  GEMS.O,  in  order  to  keep  the  generation  of 
initialseed  as  much  deterministic  as  possible,  it  has  been  chosen  to  perform  the  iterat  on 
the  individual  trajectories.  Accordingly,  to  integrate  N tj  trajectories  on  a  Grid  made  of 
Nno<fe  nodes  (Nnorfe  is  assumed  to  be  smaller  than  Ntj)  the  Master  process  (see  Fig.  1) 
generates  and  SENDS  out,  for  the  first  Nnode  iterations,  the  initial  seed  of  each  trajec¬ 
tory  without  waiting  for  the  result  of  the  integration  to  be  returned.  For  the  subsequent 
iterations  the  initialseed  of  each  trajectory  is  sent  out  only  after  one  of  the  worker  nodes 
( anynode )  has  sent  back  its  result.  After  the  sending  to  the  worker  nodes  of  initialseed  for 
the  N tj  trajectories  is  completed  the  Master  process  still  needs  to  iterate  over  the  Nno<fe 
nodes  to  collect  the  results  of  the  remaining  trajectories  and  send  a  conventional  signal 
(we  chose  to  adopt  for  this  purpose  a  negative  value  of  the  initial  seed)  to  stop  the  work 
at  worker  process  level. 

As  sketched  in  Fig.  2  the  worker  process  devoted  to  the  integration  of  each  individual 
trajectory  receives  at  the  beginning  (once  for  ever)  the  general  information.  Then  it 
iterates  on  receiving  initialseed  of  the  subset  necessary  to  generate  the  trajectory  initial 
conditions.  After  receiving  the  trajectory  initial  seed  the  worker  process,  if  the  seed  is 
not  negative,  generates  the  remaining  pseudorandom  numbers  of  the  needed  subset  to 
work  out  either  at  random  (fully  or  partially  according  to  the  chosen  distribution)  or  the 
specific  initial  value  of  the  parameters  of  the  trajectory  like  Etr,  J,  the  velocity  orientation, 
the  diatomic  internuclear  distance,  orientation  and  position  as  well  as  its  vibrotational 
(vj)  quantum  state.  The  iteration  on  time  then  starts  to  integrate  the  trajectory  from 
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Read  input  data:  v,  j,  Etr ,  A ta,  error,  maxstep,  N tj  ... 

Perform  preliminary  calculations 

BROADCAST  to  the  worker  nodes  the  needed  general  information 
Generate  initialseed  for  the  first  trajectory 

inode  0 

Do  itj  =  1,  Ntj 

inode  inode  T 1 

Generate  next  initialseed  of  the  of  the  pseudorandom  series 
IF  inode  >  Nnode  THEN 

RECEIVE  from  anynode  the  result  and  update  PVj,v'j' 
SEND  the  initialseed  to  anynode 
ELSE  SEND  the  initialseed  to  inode 

ENDIF 

EndDo  itj 

Do  inode  R  A R0de 

RECEIVE  from  anynode  the  result  and  updated  PVj,v'j' 

SEND  to  anynode  a  negative  seed  to  stop  activities 

EndDo  inode 


Figure  1:  Pseudocode  of  the  trajectory  Master  program. 


an  atom-diatom  distance  large  enough  to  consider  the  system  in  its  asymptotic  reactant 
arrangement.  The  integration  is  eventually  terminated  either  when  one  of  the  atom- 
atom  distances  has  reached  a  value  large  enough  to  consider  again  the  arrangement  as 
asymptotic  or  when  the  maximum  number  of  integration  steps  ( maxstep )  has  been  reached. 
If  at  the  integration  end  point  the  system  has  reached  the  product  asymptote  the  ,  ■/ 
probability  of  the  product  diatom  vibrotational  {v'j')  quantum  state  closest  in  energy  to 
the  computed  classical  value  is  updated. 

In  GEMS.O  the  PES  is  assumed  to  be  of  the  LEPS  type  and  the  value  of  its  parameters 
to  be  available  from  a  library  (they  may  have  been  already  calculated  using  two  other 
computational  procedures,  SUPSIM  and  FITTING,  implemented  in  our  Laboratory  (72; 
73)).  The  integration  of  a  trajectory  can  be  disposed  to  any  computing  node  of  the  grid 
while  the  integration  outcome  is  accumulated  by  the  master  process  sketched  in  Fig.  1  by 
updating  the  value  of  the  related  quasiclassical  probability  PVj,v'j'  which  can  be  displayed 
to  the  user  on  a  virtual  monitor. 

Test  runs  of  GEMS.O  (70)  performed  as  EGEE  (68)  demonstrations  have  shown  that 
the  simulator  is  highly  suitable  for  a  distribution  on  the  heterogenous  networked  envi¬ 
ronment  of  the  Grid.  In  these  runs  from  several  thousands  to  millions  or  even  billions  of 
trajectories  (to  this  end  particular  attention  has  to  be  paid  to  the  generation  of  the  pseu¬ 
dorandom  sequence)  can  be  run  depending  on  how  much  the  considered  event  is  likely 
to  occur.  The  main  feature  of  GEMS.O  is  that  of  being  cpu  bound  (this  is  in  general 
true  even  for  systems  larger  than  atom  diatom  ones).  The  implementation  of  GEMS.O 
has  impacted  EGEE  in  two  ways.  The  first  of  them  is  related  to  the  specific  requests  of 
GEMS.O  in  terms  of  infrastructure,  middleware  and  services.  The  second  is  related  to  the 
characteristics  that  molecular  simulations  need  to  possess  in  order  to  be  suitable  for  Grid 
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RECEIVE  from  the  master  node  the  needed  general  information 
[*]  RECEIVE  the  trajectory  initial seed 
IF  seed  <  0  STOP 

Generate  the  needed  subset  of  pseudorandom  numbers 
Calculate  initial  values  of  the  integration  variables 
t  =  0 

tstep  At0 

Do  it  =  1,  max  step 

t  t  -\-  tstep 

Perform  the  time  step  integration 

IF  Energy  and  total  angular  momentum  are  conserved  THEN 
IF  an  asymptotic  arrangement  has  been  reached  THEN 
perform  the  asymptotic  analysis 
SEND  results  to  Master 
GOTO  [*] 

END  IF 

ELSE  t  =  t  —  tstep  and  reduce  tstep 
END  IF 

ENDDO  it 
GOTO  H 


Figure  2:  Pseudocode  of  the  trajectory  worker  program. 


implementation. 

As  to  the  first  item  it  has  become  apparent  that  molecular  simulations  force  the  Grid 
to  better  exploit  various  levels  of  parallelization  and  distribution  using  an  appropriate 
workflow  computational  procedures  of  different  nature  and  origin.  In  other  words,  GEMS 
requests  EGEE  to  focus  more  on  an  effective  coordination  of  brainware  and  knowledge. 

As  to  the  second  item  it  has  also  become  apparent  that  EGEE  prefers  direct  calcula¬ 
tions  (with  respect  to  data  transfers)  and  coarse  grain  granularity  schemes. 
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